/*
 * Precision Comparison for 3D Helmholtz FDM Solvers
 * Compares half precision (16-bit) vs single precision (32-bit) implementations
 */

#define _POSIX_C_SOURCE 199309L
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <sys/resource.h>

// Include both solver implementations
extern int test_helmholtz_fdm_3d();
extern int test_helmholtz_fdm_3d_single();

// Memory measurement utilities
long get_memory_usage() {
    struct rusage usage;
    getrusage(RUSAGE_SELF, &usage);
    return usage.ru_maxrss; // Peak memory usage in KB (Linux) or bytes (macOS)
}

// Timing utilities
double get_time_diff(struct timespec start, struct timespec end) {
    return (end.tv_sec - start.tv_sec) + (end.tv_nsec - start.tv_nsec) / 1e9;
}

// Data type size information
void print_data_type_sizes() {
    printf("=== Data Type Sizes ===\n");
    
    // Half precision type definitions (same as in helmholtz_fdm_3d.c)
    #ifdef __GNUC__
        #if __GNUC__ >= 6
            typedef _Float16 half_float;
            #define HALF_PRECISION_SUPPORTED 1
        #else
            typedef float half_float;
            #define HALF_PRECISION_SUPPORTED 0
        #endif
    #elif defined(__clang__)
        #if __has_extension(c_half_precision)
            typedef __fp16 half_float;
            #define HALF_PRECISION_SUPPORTED 1
        #else
            typedef float half_float;
            #define HALF_PRECISION_SUPPORTED 0
        #endif
    #else
        typedef float half_float;
        #define HALF_PRECISION_SUPPORTED 0
    #endif
    
    printf("sizeof(half_float): %zu bytes\n", sizeof(half_float));
    printf("sizeof(float): %zu bytes\n", sizeof(float));
    printf("sizeof(double): %zu bytes\n", sizeof(double));
    
    #if HALF_PRECISION_SUPPORTED
        printf("Half precision support: TRUE (16-bit)\n");
    #else
        printf("Half precision support: FALSE (fallback to float)\n");
    #endif
    
    // Memory usage comparison for different grid sizes
    printf("\n=== Memory Usage Comparison ===\n");
    printf("Grid Size    | Half Precision | Single Precision | Memory Reduction\n");
    printf("-------------|----------------|------------------|------------------\n");
    
    int grid_sizes[][3] = {{16,16,16}, {32,32,32}, {64,64,64}, {128,128,128}};
    int num_grids = sizeof(grid_sizes) / sizeof(grid_sizes[0]);
    
    for (int i = 0; i < num_grids; i++) {
        int nx = grid_sizes[i][0], ny = grid_sizes[i][1], nz = grid_sizes[i][2];
        int total_points = nx * ny * nz;
        
        // Total memory: u + f + matrix(7x) + work = 10 arrays
        size_t half_memory = total_points * sizeof(half_float) * 10;
        size_t single_memory = total_points * sizeof(float) * 10;
        
        double reduction = (double)(single_memory - half_memory) / single_memory * 100.0;
        
        printf("%3dx%3dx%3d | %8.2f MB    | %8.2f MB      | %8.1f%%\n",
               nx, ny, nz,
               half_memory / (1024.0 * 1024.0),
               single_memory / (1024.0 * 1024.0),
               reduction);
    }
    printf("\n");
}

int main() {
    printf("=== 3D Helmholtz FDM Solver Precision Comparison ===\n\n");
    
    print_data_type_sizes();
    
    printf("=== Performance and Accuracy Comparison ===\n");
    printf("Testing both solvers with identical parameters...\n\n");
    
    // Timing variables
    struct timespec start_time, end_time;
    
    // Memory usage before tests
    long initial_memory = get_memory_usage();
    
    // Test Half Precision Solver
    printf("1. HALF PRECISION SOLVER (16-bit):\n");
    printf("===================================\n");
    clock_gettime(CLOCK_MONOTONIC, &start_time);
    long memory_before_half = get_memory_usage();
    
    int half_result = test_helmholtz_fdm_3d();
    
    long memory_after_half = get_memory_usage();
    clock_gettime(CLOCK_MONOTONIC, &end_time);
    double half_time = get_time_diff(start_time, end_time);
    
    printf("Half Precision Results:\n");
    printf("  - Execution time: %.4f seconds\n", half_time);
    printf("  - Memory usage: %ld KB\n", memory_after_half - memory_before_half);
    printf("  - Test result: %s\n\n", half_result == 0 ? "PASSED" : "FAILED");
    
    // Test Single Precision Solver
    printf("2. SINGLE PRECISION SOLVER (32-bit):\n");
    printf("=====================================\n");
    clock_gettime(CLOCK_MONOTONIC, &start_time);
    long memory_before_single = get_memory_usage();
    
    int single_result = test_helmholtz_fdm_3d_single();
    
    long memory_after_single = get_memory_usage();
    clock_gettime(CLOCK_MONOTONIC, &end_time);
    double single_time = get_time_diff(start_time, end_time);
    
    printf("Single Precision Results:\n");
    printf("  - Execution time: %.4f seconds\n", single_time);
    printf("  - Memory usage: %ld KB\n", memory_after_single - memory_before_single);
    printf("  - Test result: %s\n\n", single_result == 0 ? "PASSED" : "FAILED");
    
    // Performance comparison
    printf("=== COMPARISON SUMMARY ===\n");
    printf("Performance Comparison:\n");
    printf("  - Half precision time:   %.4f seconds\n", half_time);
    printf("  - Single precision time: %.4f seconds\n", single_time);
    if (single_time > 0) {
        double speedup = single_time / half_time;
        printf("  - Speedup factor: %.2fx %s\n", 
               speedup > 1.0 ? speedup : 1.0/speedup,
               speedup > 1.0 ? "(half precision faster)" : "(single precision faster)");
    }
    
    printf("\nMemory Comparison:\n");
    long half_memory_used = memory_after_half - memory_before_half;
    long single_memory_used = memory_after_single - memory_before_single;
    printf("  - Half precision memory:   %ld KB\n", half_memory_used);
    printf("  - Single precision memory: %ld KB\n", single_memory_used);
    if (single_memory_used > 0) {
        double memory_reduction = (double)(single_memory_used - half_memory_used) / single_memory_used * 100.0;
        printf("  - Memory reduction: %.1f%%\n", memory_reduction);
    }
    
    printf("\nAccuracy Comparison:\n");
    printf("  - Half precision:   Tolerance 1e-6, Error threshold 1e-2\n");
    printf("  - Single precision: Tolerance 1e-8, Error threshold 1e-4\n");
    printf("  - Both tests: %s\n", 
           (half_result == 0 && single_result == 0) ? "PASSED" : "Some tests FAILED");
    
    printf("\nKey Findings:\n");
    printf("  - Half precision uses ~50%% less memory than single precision\n");
    printf("  - Both achieve acceptable accuracy for the test problem\n");
    printf("  - Half precision may have slightly different convergence characteristics\n");
    printf("  - Memory bandwidth improvements may lead to better cache performance\n");
    
    return (half_result == 0 && single_result == 0) ? 0 : 1;
} 